Chapter 1 A New Parallel Algorithm for Computing the Singular Value Decomposition ∗

A new method is described for computing the singular value decomposition (SVD). It begins by computing the polar decomposition and then computes the spectral decomposition of the Hermitian polar factor. The method is particularly attractive for shared memory parallel computers with a relatively small number of processors, because the polar decomposition can be computed efficiently on such machines using an iterative method developed recently by the authors. This iterative polar decomposition method requires only matrix multiplication and matrix inversion kernels for its implementation and is designed for full rank matrices; thus the proposed SVD method is intended for matrices that are not too close to being rank-deficient. On the Kendall Square KSR1 virtual shared memory computer the new method is up to six times faster than a parallelized version of the LAPACK SVD routine, depending on the condition number of the matrix. 1 Background and Motivation A modern theme in matrix computations is to build algorithms out of existing kernels. In this work we propose an algorithm for computing the singular value decomposition (SVD) that requires just three building blocks: matrix multiplication, matrix inversion, and solution of the Hermitian eigenproblem. Our approach is motivated by the fact that for the Kendall Square KSR1 virtual shared memory computer there is not available, at the time of writing, a library of matrix computation software that takes full advantage of the machine. The manufacturer does, however, supply a KSRlib/BLAS library that contains a highly optimized level-3 BLAS routine xGEMM1 [6], together with an implementation of LAPACK, KSRlib/LAPACK [7], that calls the KSRlib/BLAS library. Our aim was to produce an SVD routine that is faster on the KSR1 than the KSRlib/LAPACK routine xGESVD (which implements the Golub–Reinsch algorithm), but without sacrificing numerical stability. Rather than code any of the existing parallel SVD algorithms [1], [2], [10], we wanted to make use of parallel codes that we had already developed for computing the ∗This paper is a shortened version of [4]. †Department of Mathematics, University of Manchester, Manchester, M13 9PL, England ([email protected]). The work of this author was supported by Science and Engineering Research Council grant GR/H52139, and by the EEC Esprit Basic Research Action Programme, Project 6634 (APPARC). ‡Department of Mathematics, University of Manchester, Manchester, M13 9PL, England ([email protected]). Current address: Data Information Systems PLC, 125 Thessalonikis, N. Philadelphia, Athens 142 43, Greece. This author was supported by an SERC Research Studentship. The ‘x’ in routine names stands for the Fortran data type: in this case, S or C.